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The cause of the accelerated expansion of the Universe poses one of the most fundamental ques- 
tions in physics today. In the absence of a compelling theory to explain the observations, a first task 
is to develop a robust phenomenology. If the acceleration is driven by some form of dark energy, 
then, the phenomenology is determined by the dark energy equation of state w. A major aim of on- 
going and upcoming cosmological surveys is to measure w and its time dependence at high accuracy. 
Since w^z) is not directly accessible to measurement, powerful reconstruction methods are needed to 
extract it reliably from observations. We have recently introduced a new reconstruction method for 
w(z) based on Gaussian process modeling. This method can capture nontrivial time-dependences 
in Vli{z) and, most importantly, it yields controlled and unbaised error estimates. In this paper we 
extend the method to include a diverse set of measurements: baryon acoustic oscillations, cosmic 
microwave background measurements, and supernova data. We analyze currently available datasets 
and present the resulting constraints on w{z), finding that current observations are in very good 
agreement with a cosmological constant. In addition we explore how well our method captures 
nontrivial behavior of ■w{z) by analyzing simulated data assuming high-quality observations from 
future surveys. We find that the baryon acoustic oscillation measurements by themselves already 
lead to remarkably good reconstruction results and that the combination of different high-quality 
probes allows us to reconstruct w^z) very reliably with small error bounds. 

PACS numbers: 98.80.-k, 02.50.-r 



I. INTRODUCTION 

The discovery of the accelerated expansion of the Uni- 
verse little more than a decade ago (T], [1] was a ma- 
jor surprise. Since then, many observational efforts 
to understand the underlying cause have been initi- 
ated (e.g. the Baryon Oscillation Spectroscopic Survey 
(BOSS) ll, WiggleZ % the Dark Energy Survey (DES, 
https://www.darkenergysurvey.org/), the Large Synop- 
tic Survey Telescope (LSST) Q) and proposed (e.g. Big- 
BOSS d, the Wide Field Infrared Survey Telescope 
(WFIRST) ff\, Euchd [|). These efforts focus on a set of 
diverse cosmological probes (supernovae, baryon acous- 
tic oscillations, clusters of galaxies, weak lensing, etc.) to 
combine the best possible observations in order to help 
solve this puzzle. 

The two currently most popular explanations are a 
form of dark energy or a modification of Einstein's the- 
ory of gravity on the largest observable scales. We will 
focus in this paper on dark energy as the cause for the ac- 
celerated expansion. The simplest way to realize a dark 
energy is via a cosmological constant with a dark energy 
equation of state w — p/ p = —1. A cosmological con- 
stant, however, is not theoretically well-motivated. If we 
assume that the origin is due to a vacuum energy, the 
predicted value is incorrect at the order of 10^". There- 
fore, a more natural realization of dark energy might 
be a dynamical field, similar to the inflaton that is be- 



lieved to drive the very early rapid expansion of the Uni- 
verse. Such a dynamical field, described for example by 
quintessence models [9|, would lead to a non-constant 
dark energy equation of state w(z). It is therefore one 
of the major aims of ongoing and upcoming dark energy 
missions to measure w{z) and its time variation with high 
accuracy. If w{z) is modeled via a simple parametriza- 
tion w{z) = wo — Waz/{l + z) [i3j[IHj current predictions 
for future surveys promise measurements of the constant 
part at the 1% level accuracy and of the leading time 
varying part at the 10% level. At present, the best mea- 
surements are accurate to 10% with respect to wn with 
no strong constraints on the time variation (T^ . IT3| . 

With the prospect of high-accuracy measurements 
from supernova (SN) surveys and complementary large- 
scale structure probes such as baryon acoustic oscillation 
(BAO) surveys, it is desirable to develop an accurate re- 
construction method with reliable error bars that allows 
us to extract the dark energy equation of state from dif- 
ferent measurements. While the earlier focus in the field 
was on parametric methods p^ . [Til . [l3. non-parametric 
methods are becoming more popular [l^. The major ad- 
vantage of non-parametric models is that they are not 
biased (no assumptions are made regarding the func- 
tional form for i(;(z)). A possible disadvantage might 
be that if the data quality is insufficient, non-parametric 
approaches might not provide much information about 
w{z). In principle this is also an advantage: if the data 
does not have enough information it is better to obtain 
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uncertain results with large error bars than a prediction 
which might be biased (since the functional form assumed 
for w{z) is incorrect) without this bias being reflected in 
the error bars. 

In this paper, we discuss a recently-introduced recon- 
struction method based on Gaussian process (GP) mod- 
eling [l^ [13] • A GP is a stochastic process and each 
realization is a random draw from a multivariate Nor- 
mal distribution. It is characterized by a mean and a 
covariance function, that are defined by a small number 
of parameters. Bayesian estimation methods are used to 
determine the parameters of the GP model together with 
any other physics parameters. Therefore, the final form 
of the GP model is informed by the data itself. The form 
of the covariance function is general enough to accommo- 
date a large variety of possible outcomes for wiz). The 
only assumption made is that w{z) is somewhat smooth 
and continuous. If the underlying cause for the accel- 
erated expansion is due to a physically well motivated 
reason this assumption is justified. We extend the ap- 
proach described in detail in Ref. to include differ- 
ent observational probes of u'(z), namely supernova mea- 
surements, cosmic microwave background (CMB) obser- 
vations, and BAO results. We begin with an analysis of 
currently available data. We find, not surprisingly, that 
our predictions are in good agreement with a cosmologi- 
cal constant. Using simulated data, we then explore the 
ability to extract variations of w{z) away from a cosmo- 
logical constant with improving accuracy and statistics of 
the data. The inclusion of the additional BAO and CMB 
measurements greatly help to improve these predictions. 

The paper is organized as follows. In Section |TT] we 
describe the different data sources included in our anal- 
ysis, namely, supernova, CMB, and BAO measurements. 
In Section IIIII we describe our GP model based recon- 
struction method. We carry out an analysis of currently 
available data in Section IIVI We demonstrate that the 
method will allow us to extract variations in the dark 
energy equation of state by using simulated data in Sec- 
tion |Vl Finally, we conclude in Section IVTl 



II. 



DARK ENERGY EQUATION OF STATE 
FROM DIVERSE DATA SETS 



Type la supernova measurements are currently the 
best source of information regarding possible deviations 
of w{z) from a constant value. In the future, BAO (and 
other) measurements will be a strong competitor and in 
combination will lead to the best possible constraints on 
w{z). The complementarity of the different probes is im- 
portant to break degeneracies and decrease the overall 
errors. In our previous paper [l^ we showed that our 
non-parametric reconstruction method is able to capture 
even rather sharp transitions in w{z) well if we have very 
good knowledge about flm- Supernova data alone does 
not provide this information and one needs a strong prior 
on flrn to obtain good results. This strong prior can 



be justified by the existence of complementary probes. 
In a more direct and complete implementation, multiple 
probes are included in the analysis and a joint analysis 
is performed. This allows us to relax our prior assump- 
tions on flm and to tighten our final constraints on the 
behavior of w{z). 

In the following, we provide a brief review of the differ- 
ent dark energy probes employed in this paper - super- 
novae, BAO, and CMB - and how to extract information 
about wlz) from these probes. We focus in this paper 
on the geometric probes for ■w{z). The GP analysis is in 
this case very similar for all methods - w(z) is connected 
via two derivatives with the different distance measures. 
In the next section, we explain in detail how to set up a 
joint GP model for the three different observations. 



A. Supernova Measurements 

In this paper we retain the notation from our previ- 
ous work [16 1 . For completeness, we summarize the im- 
portant equations here. The luminosity distance as 
measured by supernovae is directly connected to the ex- 
pansion history of the Universe described by the Hubble 
parameter H{z). For a spatially fiat Universe, the rela- 
tion is given by 
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(1+^)^ r 



ds 
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where c is the speed of light, TJq, the current value of the 
Hubble parameter {H{z) = a/a, where a is the scale fac- 
tor and the overdot represents a derivative with respect 
to cosmic time), and h{z) = H{z)/Hq. The assump- 
tion of spatial flatness is in effect an "inflation prior", 
although there do exist strong constraints on spatial flat- 
ness when CMB and BAO observations are combined 
(see, e.g., Ref. ^1^]). In principle, we can relax this as- 
sumption, but enforce it here to simplify the analysis. 

Instead of ^^(z), supernova data are usually specified 
in terms of the distance modulus /Lt as a function of red- 
shift. The relation between and the luminosity distance 
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10 



1 Mpc / 
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(1 + 2)C 



ds 



his) 



51ogio(i/o) + 25, 



where we used Eq. ([T]). Mb is the absolute magnitude 
of the object and ms the (_B-band) apparent magnitude. 
Writing out the expression for the reduced Hubble pa- 
rameter h{z) in Eq. ([2]) explicitly in terms of a general 
dark energy equation of state for a spatially fiat FRW 
Universe: 

h\z) ^ nr{l + Z)* + n^nil + zf (3) 

+{1 - nr - n„,){i + zf cxp (^3 rS^") 
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leads to the relation 
IJb{z) = 25 - 51ogio(7Jo) 

+ 51ogio|(l + ^)c^ ds [nr{l + sf + n,n{l + sf 



(4) 



+ (l-r!,,-r!„)(l + s)3cxp(^3^ ^ 



w{u) 



du 



-1/2 ~ 



While the term proportional to (the radiation density 
for photons and neutrinos) is negligible at low redshift, 
we include it in the equations for completeness - it will 
become important for the CMB and BAO measurements. 
We use the following relation for when CMB is 

added to the analysis: 

flriz) = r!,(0)[l + 0.227N,sf{m,a/T,)], (5) 

where for the standard three neutrino species, A^cff = 
3.04, we have TOj.a/T^ = 187/{l + z)flrh^ -10^ and /(y) ~ 
[1 + (0.3173y)^-*^]^/^-^^ (In the following, wherever we 
quote values for fir they are quoted at z = 0.) 

Note that Hq in Eq. (U) cannot be determined from 
supernova measurements in the absence of an indepen- 
dent distance measurement. Thus Hq can be treated as 
an unknown and absorbed in a re-definition of the ab- 
solute magnitude Mb = Mb — Slog^giJo + 25, which 
accounts for the combined uncertainty in the absolute 
calibration of the supernova data, as well as in Hq. Us- 
ing this, the _B-band magnitude can be expressed as 
niB — 5 logiQ Dl{z)+M b where Dl{z) = HodL{z) is the 
"Hubble-constant-free" luminosity distance (throughout 
this paper we will follow the convention to use capital 
letters for Hubble-constant-free distances and small let- 
ters for distances measured in Mpc. Different papers use 
different conventions.). The measurement of /ub is only 
a relative measurement and M b allows for an additive 
uncertainty which can be left as a nuisance parameter. 
To simplify our notation, we absorb 5 log]^Q(iJo) ~ 25 into 
our definition of the distance modulus, leading to: 

MB =MB + 51ogio(i^o)-25-51ogio[i5L(;2)]. (6) 

With this definition of the distance modulus we have cal- 
ibrated the overall offset of the data to be zero. To ac- 
count for uncertainties in this calibration, we introduce 
a shift parameter with a broad uniform prior. The 
expected value for = 0. 



B. BAO Measurements 

Baryon acoustic oscillations provide another powerful 
measurement of the expansion history of the Universe. 
Similar to supernovae they yield a geometric probe of 
dark energy. By carrying out measurements of the clus- 
tering along the transverse BAO scale one can obtain the 
angular diameter distance d^(z), defined as 



d-Aiz) 



1 



ds 



1 + zHoJo h{s) 



(7) 



and by measuring the BAO scale along the line of 
sight, one obtains information on the Hubble parame- 
ter H{z) itself (for details on future measurements, see, 
e.g. Ref. [gI]). Both of these measurements will be carried 
out in terms of the sound horizon at the epoch of baryon 
drag, rs{zd), given by: 



ds 



H{s)Jl 



3^6 



(8) 



and the final measurements will be in terms of dA{z)/rs 
and H(z)rs. Current data provide information only on 
the angular diameter distance. The structure of Eq. ([7]) 
with respect to its w-dependence via two integrals is ex- 
actly the same as for (ii(z) given in Eq. ([T]). This makes 
it very easy to carry through a reconstruction approach 
combining both probes. 



C. CMB Measurements 

For the CMB measurements we employ the so-called 
shift parameter i?(z*) first introduced by Bond et al. jli| : 



R{z 



*) = 1 + zAdAiz^, 



h{sy 



(9) 



where is the redshift of decoupling (z* ~ 1090) and 
the angular diameter distance dA is given in Eq. ([7]) . The 
shift parameter is related to the peak heights and the lo- 
cations of the peaks in the temperature power spectrum 
of the CMB. As we will show in our analysis below, the 
shift parameter is very helpful in breaking the degeneracy 
between and w{z) when used in a combined analysis 
with supernova data. As an alternative to using the full 
CMB power spectra, the shift parameter provides a good 
way to summarize (see, e.g., Ref. [l^]) CMB measure- 
ments, hence simplifying dark energy investigations. 

One caveat of using R as pointed out in, e.g., Ref. [18 1 
is the fact that i? is a derived quantity from fitting to the 
CMB power spectrum and therefore assumes a certain 
cosmology. It is therefore important to state explicitly 
the assumption made under which the best-fit value for 
R was derived. Several groups including Refs. (20l - [2^ 
have studied this point in more detail and found that 
the constraints on R are relatively stable under minor 
modifications of the dark energy parameters underlying 
the analysis, including dark energy clustering (22| . It 
was found that massive neutrinos had a larger effect on 
R (few percent level) [l^]. In Ref. ^23] an analysis of 
WMAP-3 data was carried out and it was found that for 
non-flat cosmologies, the value for R was very similar for 
different dark energy models, including constant w and 
time- varying w parametrized via wq + Wa{l — a). In ad- 
dition, the best-fit values for R in the current WMAP-7 
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analysis are the same within error bars for different un- 
derlying cosmologies, including wCDM and open wCDM 
models. 

In our analysis of currently available data it should 
be kept in mind that we use the best-fit value for R de- 
rived under the assumptions of a flat FRW universe with 
w = — 1, an effective number of neutrinos of N^ff = 3.04 
and a primordial power spectrum close to a power law. 
As we show below, the inclusion of R in the analysis 
in addition to the supernova data does not alter the re- 
sult for w{z) itself, its main contribution is to help relax 
the assumption on $7^- For this reason, the fact that 
the value we use for R is derived for a specific model is 
of much less consequence. In the case of our simulated 
data, the value of R is obtained for the correct underly- 
ing cosmology, in which case the above discussion does 
not apply. 

Another issue arises with the CMB measurement point 
due to its origin at high redshift. The SNe and BAO data 
points occupy a redshift range between z £ (0, 2) making 
it easier to set up a coherent non-parametric reconstruc- 
tion approach. The CMB data point on the other hand 
is a single point around z ~ 1000, so far away that it is 
bound to cause problems for any non-parametric method. 
Consequently, we have to make some assumptions about 
the behavior of w{z) in the range z G (2,oo) - the sim- 
plest choice is w = const. 



III. RECONSTRUCTION WITH GAUSSIAN 
PROCESS MODELING 

A. Overview 

We have recently introduced a nonparametric recon- 
struction method based on GP modeling and Markov 
chain Monte Carlo (MCMC), and applied it to super- 
nova data [l^, [l3] . We refer the reader to these papers 
for details on the implementation of the GP model. Here, 
we provide a general introduction and explanation of the 
idea behind the reconstruction process with GP models 
and then focus on how to extend the method to include 
multiple data sources. 

Gaussian processes extend the multivariate Gaussian 
distribution to function spaces, with inference taking 
place in the space of functions. The defining property of 
a GP is that the vector that corresponds to the process at 
any finite collection of points follows a multivariate Gaus- 
sian distribution. Gaussian processes are elements of an 
infinite dimensional space, and can be used as the ba- 
sis for a nonparametric reconstruction method. Gaussian 
processes are characterized by mean and covariance func- 
tions, defined by a small number of hyperparameters [23j . 
The covariance function controls aspects such as rough- 
ness of the candidate functions and the length scales on 
which they can change, aside from this, their shapes are 
arbitrary. The use of Bayesian estimation methods (in- 
cluding the MCMC algorithm) allows us to estimate the 



hyperparameters of the GP correlation function together 
with any other parameters, comprehensively propagat- 
ing all estimation uncertainties [2J| . Using the definition 
of a GP, we assume that, for any collection zi,...,Zn, 
w{zi), ...,'w{zn) follow a multivariate Gaussian distribu- 
tion with a constant negative mean and exponential co- 
variance function written as 



K{z,z') = 



(10) 



Here p G (0, 1) is a free parameter that, together with k 
and the parameters defining the likelihood, are fit from 
the data (p and k are the hyperparameters of the GP 
model). The form of the assumed correlation function 
implies that, theoretically, there is non-zero correlation 
between any two points. The parameter p controls the 
exponential decay of the correlation as a function of dis- 
tance in redshift, but it does not provide a bound for the 
correlation between two points. 

The value of a G (0, 2] influences the smoothness of the 
GP realizations: for a = 2, the realizations are smooth 
with infinitely many derivatives, while a = 1 leads to 
rougher realizations suited to modeling continuous non- 
differentiable functions. Here we use a = 1 to allow for 
maximum flexibility in reconstructing w. (For a com- 
prehensive discussion of different choices for covariance 
functions and their properties, see Ref. [23|].) We set up 
the following GP for w. 



w{u) - GP {i),K{u,u' j). 



(11) 



The process is started using a mean value of •& = — 1; 
given current observational constraints on w this is a 
natural choice. Even though the mean is fixed, each GP 
realization actually has a different mean with a spread 
controlled by k and the means are adjusted during the 
analysis to slightly different values suggested by prelimi- 
nary runs (we use this strategy for some of the simulated 
data sets below). This adjustment is purely informed by 
the data and demonstrates the flexibility of the approach. 
In principle, the mean could also be left as a free parame- 
ter. After the adjustment we measure the posterior mean 
and ensure that it is close to the prior mean. 



B. Combining Multiple Data Sources 

In order to determine the optimal values for the GP 
modeling parameters and the cosmological parameters, 
we follow a Bayesian analysis approach .25] . We use 



MCMC algorithms to fit for the parameters [2j|, result- 
ing in posterior estimates and probability intervals for 
^rn and A^, and the hyperparameters that specify the 
GP model, k and p. We choose the following priors for 
the hyperparameters: 



7r(«)^/G(6,2), 
7r(p) - Beta{6, 1). 



(12) 
(13) 
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Here the notation "~" means "distributed accord- 
ing to" , and IG is an inverse Gamma distribution 
prior, with the probabihty density function f{x;a,f3) — 
/3"2;-"-ir(a)-i exp(-/3/x), with x > 0. The probabil- 
ity distribution of the Beta prior is given by f{x; a, /3) = 
T{a+P)x"-\l-xy-^/\r{a)T{(3)] (for examples of these 
distributions, see, e.g. [17]). 

Turning to the cosmological parameters, we choose: 

7r(r2„) - 7V(0.27, 0.04^) SN data only, (14) 

7r(rim) ~ U{0, 1) combined analyses, (15) 

7r(A^)^C/(-0.5,0.5), (16) 

7r(cr2) oc a-^ SN data, (17) 

7r(cr|) cx BAG data, (18) 

where C/ is a uniform prior, with the probability den- 
sity function f(x]a,h) = \/{b ~ a) for x G [a, 6] and 
otherwise. is a Gaussian (or Normal distributed) 
prior with the probability density function /(a;;/i,CT^) — 
exp[— (a; — /^)^/(2t7^)]/v27R7^. The squared notation for 
the second parameter in N{^, ct^) is used to indicate that 
a is the standard deviation (to prevent possible confusion 
with the variance a^). (The parameters in the U and IG 
distributions do not have this same meaning of mean and 
standard deviation as in the Normal distribution.) 

The prior for VLm for the analysis of supernova data 
alone is informed by the 7- year WMAP analysis [3l for 
a wCDM model combining CMB, BAG, and Hq mea- 
surement. Gnce a second cosmological probe is included 
in the analysis, the assumption on this prior can be re- 
laxed and we choose a uniform prior for the analysis of 
the combined data sets. We also allow for an uncertainty 
in the overall calibration of the supernova data, A^, and 
choose a wide, uniform prior for A^. 

Next we discuss the likelihoods for the different probes. 
We assume that the SNe, CMB, and BAG measurements 
are independent of each other which allows us to derive 
a likelihood for each probe separately. The likelihood for 
the supernova data is given by: 



correlated and bivariate Normal distribution, given by: 



L=^^(a,0)oc 



exp 



2 ^ 

i=l 



(19) 

where 9 encapsulates the cosmological parameters as well 
as the hyperparameters, i.e., {A^, 17^, p\ and cr^ is the 
associated variance, expected to be close to unity. For the 
CMB data we have an equivalent expression: 



1 



l(y*-R{z*, 




(20) 



Since we only have one data point, we cannot assign a 
variance parameter. The likelihood for the BAG data is 
slightly more complicated. For each BAG point we have 
two observed distance measures. These measurements 
{yiii y2i) are correlated and we assume that they have a 



where 



y2t 



K = 



MVN 



DA{zi)/rs 
H{zi) * rs 



,a|K 



(21) 



(22) 



and gb is the associated variance parameter. This leads 
to the following likelihood for the BAG data: 



BAO 



L 



with 



{aB,0) (X 



1 



(7|K|™/2 



exp 



y2i - H{z) * Ts 



(23) 



(24) 



We can find the combined likelihood simply by multiply- 
ing the likelihoods since we assume the different probes 
are uncorrelated: 



jtot&\ _ 2^SNo ^ ^CMB ^ ^BAO 



(25) 



IV. RESULTS FOR CURRENT OBSERVATIONS 

We begin our analysis by reconstructing w{z) from cur- 
rently available data. We use the supernova data set re- 
cently released by Amanullah et al. [ij]. This so-called 
Union-2 compilation (extending the Union compilation 
from Ref. u^) consists of 557 supernovae between red- 
shift z = 0.015 and z = 1.4. The magnitude errors in the 
data set range between 0.08 and 1.02, with an average 
error of ^ 0.2. 

In addition to the supernova data we include the most 
recent BAG measurements from the Two-degree-Field 
Galaxy Redshift Survey (2dFGRS}^ at z = 0.2 and the 
Sloan Digital Sky Survey (SDSS) j26j] at z = 0.35 given 

by 

rs{zd)/dv{z = 0.2) ^ 0.1905 ± 0.0061, (26) 
rs{zd)/dv{z = 0.35) 0.1097 ± 0.0036, (27) 

where dv{z) = [(1 + zfd\cz/H{z)Y/^. For the CMB 
analysis, we use the most recent measurement of the shift 
parameter R from WMAP- 7 [l^, given by 



R{z^) = 1.719 ±0.019. 



(28) 



In order to have a complete description of the problem we 
have to specify some additional cosmological parameters 
that are expected to have little or no effect on dark en- 
ergy. These parameters - fixed at the best-fit WMAP- 7 
values from their ACDM analysis - are: $1^ = 4.897-10^^, 
Zd = 1020.3, z^ = 1090.79, and rj^/fl^ = 914.54. 



6 



SNe 


SNe+CMB 


SNe+BAO 


SNe+BAO+CMB 










0.2 0.4 0.6 0.8 1.0 1.2 0.2 0.4 0.6 0.8 1.0 1.2 0.2 0.4 0.6 0.8 1.0 1.2 0.2 0.4 0.6 0.8 1.0 1.2 
z z z z 




0.0 0.2 



0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0 

^Vn ^ ^Kn 



FIG. 1. Results for reconstructing ■w{z) from currently available data. The top row shows reconstruction results for w{z) (red 
line; the black dashed line shows w — —1) for an exponential covariance function, i.e. q = 1, including different cosmological 
probes, the second row shows the corresponding posterior for Jl^ (red lines: priors, black lines: posteriors). The first column 
shows the results from supernova data only, the second column includes CMB measurements, the third column uses supernova 
and BAO data, and the fourth column shows the results for a combined supernova-BAO-CMB analysis. The light blue contours 
show the 95% confidence level, the dark blue contours the 68% confidence level. As is to be expected, the error bars shrink 
somewhat if more data sources are included, though the effect is small due to the limited number of extra data points. As 
found previously by others (e.g. [13]), current data are consistent with a cosmological constant. 



TABLE I. Union 2 Data set- 95% Pis, the last two columns are results from Ref. [13], Table 11 for comparison. 
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We carry out four different analyses: supernova data 
by themselves with a Gaussian prior for f2,„ given in 
Eq. (|14p . and combined analyses for supernova data and 
CMB, supernova data and BAO, and for all three probes. 
For the combined data sets we can relax the prior as- 
sumptions on flm and use a wide uniform prior, given 
in Eq. (llSp . The results are summarized in Table |T] and 
Fig.ffl 

All results are consistent with a cosmological constant, 
i.e. w = —1, as can be seen in the first row in Fig. [T] 
The supernova data by themselves have no constraining 
power on f2m and therefore force us to choose a rather 
strong prior. The lower panels in Fig. [T] show the prior 
(red line) and posterior (black line) for 51^ demonstrating 
this point clearly. If we include either CMB or BAO or 



both, the constraints on 51„i get much better. As can 
been seen in Table HI the error estimates for Qm shrink by 
almost a factor of two if all probes are combined. Overall, 
the supernova data by themselves lead to a slightly higher 
value of rim , while the combination with CMB data leads 
to a lower value. The inclusion of the BAO points shifts 
up the value for r2„j considerably, by more than 10% 
compared to the supernova-CMB analysis. Nevertheless, 
within the error bars, all values for flm are consistent and 
agree well with the best-fit WMAP-7 values including 
different probes. The value for the shift parameter is 
very close to zero in all cases. 

A brief comparison with Ref. jT3| also shows very good 
agreement. For ease of comparison, we quote their results 
in the last two columns of Table |T] for the case of a flat 
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Universe and w = const. The trends in the best-fit value 
for are exactly the same as we find, the value is low- 
est for the case of supernova-l-CMB data and highest for 
supernova-|-BAO data. They also find that for the super- 
nova -fBAO analysis, w is slightly below w = —1 while 
for all other cases it is very close to w = —1. Their anal- 
ysis is also consistent with a cosmological constant. It 
is very interesting to note that our error estimates for 
w{z) are also very similar to the findings of AmanuUah 
et al., even though their assumption ofw — const, is very 
restrictive. This is very encouraging since it shows that 
our method leads to tight error bounds without loss of 
flexibility in allowing for time variations in w{z). In con- 
trast, a Wo ~ Wa fit would have increased the error bars 
considerably. 



model with w{z) = -1.0006 -I- 308472/(exp[20/(l -t- z) -I- 
617439]). The resulting equations of state are shown in 
Figure [21 For each model we choose ftm — 0.27 and fix 
Ho = 70.4 km/s/Mpc, 0.0226, ujr = 2.469 • 10-^ 

= 1090.89, and za = 1020.5. While Model 3 is already 
ruled out observationally it provides a good example for a 
rather sharp transition in w{z). For each model we create 
two data sets: (i) We assume the best-possible scenario, 
a space mission to obtain supernova measurements out 
to redshift z = 1.7 and in addition a BigBOSS-like BAO 
survey, and CMB data; (ii) good ground based supernova 
measurements in combination with BigBOSS and CMB 
measurements. In the following we provide some details 
on the assumptions for the different data sets. 



V. RESULTS FOR SIMULATED DATA 

In this section we investigate how well our method 
works for reconstructing w(z) with future high-quality 
data. Current limitiations - uncertainties in the data 
and limited statistics - prevent us from extracting possi- 
ble time variations in w(z) reliably. The error bands are 
still rather large and results are in complete agreement 
with a cosmological constant. Future measurements will 
hopefully change this: if there is a small time variation 
in w{z) we should be able to detect it. In our previous 
paper [l6| we generated a supernova data set, assuming 
high-quality measurements from a WFIRST-like mission. 
We showed that a set of ^ 2300 supernovae out to a red- 
shift of z — 1.7 and perfect knowledge of flm allows us 
to confidently extract time variations in the dark energy 
equation of state. We also showed that larger uncertain- 
ties in flm degraded this result due to well-known degen- 
eracies between w and flm- These degeneracies can be 
broken by including different data sources. We show in 
the following that a combination of accurate supernova 
data with results from a BAO survey such as BigBOSS 
will provide sufficient information to enable a reliable 
and interesting reconstruction of 'w{z). The combination 
of these different data sources eliminates the degeneracy 
problem and provides reliable constraints on the time 
variation of w{z) without requiring "perfect" knowledge 

of ^Im- 



A. The Simulated Data 

We generate simulated data for all three probes (su- 
pernovae, BAO, CMB) and three different cosmological 
models. The models used are the same as in our pre- 
vious work (see Ref. [l6i for more details). Model 1 
has a constant dark energy equation of state w — —1, 
Model 2 is based on a quintessence model with a min- 
imally coupled scalar field and a dark energy equation 
of state w{z) = (0/2 - Vb<^2)/(<^72 + Vocj}"^), and for 
Model 3 we choose a slightly more extreme quintessence 



1. Supernova Measurements 

As mentioned above we investigate two different sets 
of simulated supernova measurements. The first one is 
the same as we used in Ref. 0|. It contains 2298 data 
points distributed over a redshift range of < z < 1.7 
with larger concentration of supernovae in the midrange 
redshift bins (0.4 < z < 1.1) and at low redshift (z < 
0.1). The exact distribution is shown in Ref. [l^ in Fig.l. 
For the distance modulus we assume an error of = 0.13. 
The measurements are presented in the following form: 

= ai + e,. (29) 

In this notation, the observations jli follow a normal dis- 
tribution with mean a{zi), the standard deviation being 
set by the distribution of the error, e^, representing a 
mean-zero normal distribution with standard deviation, 
TiU. Here, is the observed error and a accounts for a 
possible rescaling. In addition, we assume that the errors 
are independent. 

For the second set of simulated supernova data we con- 
sider the same number of data points as currently avail- 
able from ground-based surveys (557 measurements). 
The redshift distribution is shown in Fig. [3] The dis- 
tribution extends to z = 1.4 with a maximum at low 
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FIG. 2. Dark energy equation of state w[z) for our three 
simulated models. 
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FIG. 3. Redshift distribution of the small supernova data set. 
The simulated data has exactly the same distribution as the 
real data. 
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FIG. 5. Two realizations of 20 simulated BAO points for a 
BigBOSS-like survey. Left column: angular distance diam- 
eter, right column: Hubble parameter both in terms of the 
sound horizon. Model 1 {w — const.) is shown with error bars 
with one standard deviation. For Model 2 (green dashed) and 
Model 3 (blue dotted) we show the exact predictions. 



3. BAO Measurements 



FIG. 4. Small supernova data set with error bars and the 3 
simulated models. 



redshift and around z = 0.3. Only a handful of super- 
novae are available at higher redshifts. Since we assume 
that the measurements are taken from the ground, we 
increase the errors on the distance modulus to = 0.15. 
Figure |4] shows the distance modulus redshift relation of 
these measurements for Model 1 with error bars. The 
exact relations for Model 2 and 3 are shown in addition. 
The differences between the three models are very small, 
pointing to the challenge of the reconstruction task. 



2. CMB Measurements 

For the CMB points we use the following realizations 
(the exact values for R for each model are given in paren- 
theses): 

Model 1 :i?(20 = 1-736 ± O.G19(i?°'' = 1.723), (30) 
Model 2 :i?(z*) = 1.716 ± 0.019(i?'='' = 1.702), (31) 
Model 3 :i?(z*) = 1.683 ± 0.019(i?'='' = 1.670). (32) 



Future BAO surveys such as BigBOSS will obtain 
measurements of the angular diameter distance, dA{z), 
as well as the Hubble parameter H{z), in terms 
of the sound horizon at the epoch of baryon drag, 
rs{zd)- For our simulated BAO data sets, we follow 
the specifications for a BigBOSS survey as outlined in 
http://bigboss.lbl.go v/docs/BigBOSS-NOAO -public.pdfl 
We assume a survey area of 24000 deg^ (covering north- 
ern and southern skies) and adopt the galaxy density 
distribution estimated in the BigBOSS proposal (Table 
2.3 in the aforementioned document). In this proposal, 
measurements from luminous red galaxies and emission- 
line galaxies are combined. The resulting distribution 
accounts for several sources of inefficiency (discussed 
in the BigBOSS proposal) leading to a degradation of 
the galaxy number density at high redshift. Often, a 
constant galaxy density over the whole redshift range is 
assumed. We studied this case as well and found that 
the eventual results in both cases are very similar. In 
order to derive estimates for the errors of the simulated 
measurements, we use a publicly available code intro- 
duced in Ref. [2^. The formula used to obtain BAO 
errors in this code is a 2D approximation of the full 
Fisher matrix formalism. In j28|, the results for the full 
Fisher matrix calculation and this method are shown to 
match well. Although these results are for ACDM, they 
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should hold for other cosmologies. 

The input parameters for the code are: cts at the 
present epoch, 'Ej_ = T,qG = transverse rms Lagrangian 
displacement, with G — growth factor normalized such 
that G = (1 + z)-^ at high redshift, Eg = UAh'^ Mpc 
for a cosmology with erg — 0.9 at present and scaling 
linearly with erg; S|| — T,qG{1 + /) = line of sight rms 
Lagrangian displacement, with / — (i(lnG)/(i(lna), G, Sq 
as before; and the number density = 3 x 10~^ft.^/Mpc^ 
(0). G, /, CTs are input correctly for each model. The 
biggest possible source of error are the formulae used for 
I]_L, ; these were shown to be reasonable fits to the true 
values in [2^. The value of Eg given is also for the cos- 
mology used in [2^ . For a different cosmology, Eq would 
obviously be different, and the simplest way to deal with 
this, as suggested in the paper, is to scale it linearly with 
as- This may not be completely accurate as we use very 
different cosmological models but should yield a reliable 
estimate. 

Figure [5] shows two realizations for a ACDM model 
(Model 1) for the angular distance diameter Dji{z)/rs in 
the left column and for the Hubble parameter H{z)rs in 
the right column. In addition, we show the exact pre- 
dictions for Model 2 and 3. We use two realizations for 
the BAO data to demonstrate the dependence of the re- 
construction as a function of realization. Because obser- 
vations represent only one realization, this imposes an 
an irreducible limitation on the reconstruction program, 
whether non-parametric or not. We will return to this 
issue in future work. 



B. Results 

1. Prelude 

Before we present our results for the combined analysis 
of different cosmological probes we show the constraints 
we obtain from the simulated BAO data alone on w{z). 
The results are already remarkably good. We choose a 
flat prior for VL„i for this analysis. 

Figure [B] shows the results for both realizations pre- 
sented in Figure [SJ Table |lT] provides the best fit values 
for Vim all three models for the left column (Realization 
I) and Table Hill for the right column (Realization II). For 
Model 1 (first row) the predictions are slightly low for the 
first realization but overall the results are consistent with 
the input model., w — We verified that this result 
does not change considerably if the tighten the prior on 
Vim- Similar trends can be seen for Model 2 and 3. We 
will come back to these trends later in the discussion on 
the results for combined data sets. The value for 
for realization I (Table |ll]) is slightly high in all cases - 
adding CMB measurements decreases the error on Vim 
but in fact shifts the best fit values even higher. The 
second realization leads to values for Vim very close to 
the input value for BAO measurements only, the CMB 
point again shifts it up slightly. The reconstruction from 



the BAO data only works remarkably well - in all cases 
the underlying model is captured within the error bars 
reliably. 



2. Combining Different Data Sets 

Next we present the results for Model 1-3 for several 
different combinations of data as discussed above: 

Ground-based supernova mission (Figs. upper rows; 
Table\n\): 

• 557 supernovae out to z = 1.4, = 0.15 
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FIG. 6. Left column: reconstruction result for wi^z) for re- 
alization I of the BAO data later used in combination with 
the small supernova sample. Right column: results for real- 
ization II later used in combination with the large supernova 
sample. Top to bottom: results for Models 1-3. The dashed 
line shows the underlying theoretical model, the dark blue 
region shows the 68% confidence level, the light blue region 
the 95% confidence level, the dark blue line shows the mean 
reconstructed history. In all cases, the reconstruction results 
capture the "truth" within the error bands reliably. 
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FIG. 7. Reconstruction results for the model with w = —I. From left to right different probes are considered, SNe, SNe+CMB, 
SNe+BAO, and a combination of all three measurements. The red dashed line shows the truth, the blue solid line the mean 
result for the reconstruction. The blue shaded region shows the 68% confidence level, the light blue shaded region the 95% 
confidence level. The upper row shows the results for the small supernova data set (557 supernovae with = 0.15 out to 
z — 1.4) while the lower row shows potential space-based supernova measurements (2298 supernovae with = 0.13 out to 
z — 1.7). The CMB data point is the same in all cases where it is included, the BAO data are of same quality but two different 
realizations out to 2 = 2. Note that the redshift range varies in the different panels depending on which probes are included. 



• supernovae + CMB measurement 

• supernovae + 20 BAO points (realization I) 

• supernova + BAO + CMB measurements 

Space-based supernova mission (Figs. 00 lower rows; 
TableUEl: 



• 2298 supernovae out to z = 1.7, t; — 0.13 

• supernovae + CMB measurement 

• supernovae + 20 BAO points (realization II) 

• supernova + BAO + CMB measurements 

As for the real data, we choose a stronger prior for 
flm in the case of analyzing supernova data only while 
we use a flat prior for any combination of data. Fig- 
ure [7] shows the results for Model 1. The reconstruction 
from supernova data only works very well - the additional 
data points (comparing the upper and lower panel) help 
reduce the error bands (note that the redshift range in 
the lower row showing the results for 2298 supernovae 
extends out further) and also lead to a better estimate 



for with tighter error bounds, given in Tables HIl IIIII 
The addition of the CMB point (second column in Fig- 
ure [7]) allows us to choose a much less strict prior on firm 
i.e. a flat prior. Overall, the reconstruction works well 
with the combination of supernova and CMB measure- 
ments, the error bands on w{z) shrink considerably. The 
estimate for ri„i is slightly too high leading to a small 
overall underestimation of w{z) [we remind the reader of 
the degeneracy of flm and w{z)]. In the third column we 
show the supernova+BAO analysis. In this case, both 
results extend to z = 1.7 due to the BAO data at those 
redshifts. In the upper row, the supernova data only cov- 
ers a redshift range out to z = 1.4, the overall result is 
similar to the result from the BAO data only (Figure [5]) 
though the error bands shrink considerably. Combining 
all three data sets leads to even narrower error bands 
(fourth column). In the lower row the small downward 
trend from the CMB point is compensated by the small 
upward trend from the BAO measurements at high red- 
shifts, leading to an almost perfect reconstruction result. 
In the upper row, both CMB and BAO realization have a 
small downward trend in w{z) which surveys in the final 
result. Overall, the "truth" is captured well in all cases 
and lies well within the error bounds. We would like to 
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FIG. 8. Same as in Fig. [7] but for Model 2, the quintessence model. 



TABLE II. Posterior 95% Pis, 557 SNe points 
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emphasize that the dark blue line in the figures only rep- 
resents the mean of the reconstruction result; much more 
significant are the error bands themselves - these must 
capture the true underlying model to establish a valid 
approach. 

The results for Model 2 and 3 are similar, shown in 
Figures [5] and ini Model 2 exhibits a small time variation 
which could be extracted from future data. The powerful 



combination of all three probes can be gauged by the 
relatively small error bands shown in the fourth column 
in Figure [51 At low to intermediate redshifts (out to 
z ~ 0.6) a cosmological constant is clearly disfavored. 
The supernova data alone would not have had enough 
information to disfavor w = — 1 at any redshift, as the 
error bands in this case clearly include a cosmological 
constant. The inclusion of high redshift supernova data 
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FIG. 9. Same as in Fig. [T] but for Model 3. 



TABLE III. Posterior 95% Pis, 2298 supernova 
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improves the results somewhat, the overall reconstruction 
shown in the lower left corner of Figure |8] is excellent with 
narrow error bands. In this case, the constraints for f2„ 
are also very close to the input value for the theoretical 
model with tight error bands. 

Model 3 has a rather strong variation in w{z). While 
this model is observationally ruled out already, it pro- 
vides a good test bed for our new approach to demon- 



strate that more complicated dark energy equation of 
states can be reconstructed. As we discussed in detail 
in Ref. [l^ the degeneracy between fi^ and w{z) makes 
the reconstruction task rather difficult - the left panels in 
Figure |9] show the constructed w{z) from supernova data 
only with a Gaussian prior on VL^n- The error bars are 
rather wide and include a cosmological constant comfort- 
ably. The addition of the CMB point already improves 
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the result considerably, in this case we choose a flat prior 
on Ujn- The best-fit value for is very close to the 
input value of 0.27 compared to the case where we ana- 
lyze supernova data only. The inclusion of the BAO data 
(third and fourth column) in both cases (557 and 2298 
supernova data points) improves the results even more. 
The time dependence is well captured and the estimate 
for rtm is also very good. 

Some final remarks on the content of Tables HIl IIIII in 
addition to the results discussed above, we provide some 
information on the results for the combination of BAO 
and CMB measurements. Overall, the extra information 
from the CMB measurement does not help very much to 
improve the results, contrary to what we find when we 
add this information to the supernova data. In addition 
to the constraints on the cosmological parameters and 
error behavior of the data (given by a for the supernova 
data and (Jb for the BAO data) we list the final hyper- 
parameters for the GP model in the last three columns. 
Perhaps the most interesting parameter here is the ad- 
justed mean value for w(z) given by in the last column. 
As we described in Ref. [la| in detail, we start the GP 
model with some value for i) (in the case of Model 1, 
■i? = — 1 is the natural choice for example) and run the 
reconstruction program for some time. The results then 
have information about an improved value for the mean 
of the GP model and the analysis framework can be ad- 
justed accordingly. As can be seen in the Tables, the final 
values for d are close to the mean value of the underlying 
truth. Because the adjustment scheme works extremely 
well, we started basically all reconstruction evaluations 
at 1? = — 1, the GP model automatically suggesting bet- 
ter mean values if the choice was non-optimal. Overall, 
the reconstruction of w(z) works very well when multiple 
sources are included. 



VI. CONCLUSION 

In this paper we have introduced a new non-parametric 
reconstruction scheme for the dark energy equation of 
state w(z) combining multiple cosmological probes. The 
reconstruction scheme is based on a GP modeling ap- 
proach and provides very good constraints on w{z) with 
reliable error bars. The basic method was introduced in 
Ref. [ll] for supernova data only. Here we extend the 
methodology to include BAO and CMB measurements. 



We have carried out an analysis of currently available 
data and found excellent agreement with a cosmological 
constant consistent with a large number of recent pub- 
lications, including Refs. [H, [ll, HI, [H, [l^]. We have 
also demonstrated our method on simulated data for dif- 
ferent cosmological models. In all cases, the GP model 
approach performed very well. 

An important aspect of our new approach (as stressed 
in Refs. [l^ [l3l) is the simultaneous constraint of the 
cosmological parameters as well as the hyperparameters 
of the GP model from the data. In comparison to para- 
metric approaches, our new method is more flexible and 
can therefore capture even subtle time variations in w{z) 
if the data quality is good enough. It produces narrow 
error bands over the full redshift ranges considered. For 
a more detailed comparison with parametrized methods, 
see Ref. [11]. 

The combination of different data probes mitigates the 
problem of degeneracies between w{z) and flm as to be 
expected. An encouraging observation is that even the 
BAO data alone (of high quality from a BigBOSS-like 
survey) can deliver good constraints on the time depen- 
dence of the dark energy equation of state, clearly com- 
petitive with space based supernova observations. 

Our new non-parametric reconstruction approach 
lends itself to analysis of the promise of future dark en- 
ergy probes in a reliable way. For example, possible ten- 
sion in the data due to e.g. insufficient understanding of 
systematic errors would lead to an increase in the error 
bands when combining different probes (a different at- 
tempt to solve this problem with parametric methods is 
discussed in e.g. Ref. 0). The GP based approach can 
therefore help to optimize future dark energy missions. 
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